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We developed a genetic algorithm (GA) in the Heisenberg model that combines a triadic 
crossover and a parameter-free genetic algorithm. Using the algorithm, we examined the 
ground-state stiffness of the ± J Heisenberg model in three dimensions up to a moderate size 
range. Results showed the stiffness constant of = in the periodic-antiperiodic boundary 
condition method and that oi 9 ^ 0.62 in the open-boundary-twist method. We considered 
the origin of the difference in 9 between the two methods and suggested that both results show 
the same thing: the ground state of the open system is stable against a weak perturbation. 
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1. Introduction 

Optimization methods have found widespread application in computational physics. 
Among these, the investigation of the low-temperature behavior of spin glasses (SGs) has 
attracted much attention within the statistical physics community because, despite its simple 
definition, its behavior is far from understood. Particularly, domain-wall energies of the SGs 
at the absolute zero temperature (T = 0) are suggested to give a stable SG phase at a very 
low temperature. From a computational point of view, the calculation of spin-glass ground 
states is very demanding and various algorithms have been developed.^ Pal found that ge- 
netic algorithms (GAs) are useful in searching for the ground state of the ±J Ising model. ^ 
He proposed a triadic crossover together with an effective local optimization method. The 
GA algorithm has been improved in several ways; a powerful local optimization method is 
attached,^ and the GA is coded in a different manner.^' ^ The system size that is treatable 
in a 3D SG model is moderate at 13 x 13 x 13. Unfortunately, those algorithms more or less 
use the nature of the Ising spin and do not apply continuous spin models. One of the authors 
and his coworkers showed that a GA also works in the Heisenberg SG model.® They applied a 
one-point crossover together with interface optimization. Although it enables us to treat the 
3D lb J model up to a moderate size, its coding was rather complicated; several improvements 
are necessary for application to different systems. 

In this paper, we will develop a GA in a different manner for the Heisenberg model and 
consider the ground-state properties of the 3D it J Heisenberg model on L x L x L described 
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by the Hamiltonian 

H = - ^ JijSi ■ Sj, (1) 

(ij) 

where Si is the Heisenberg spin of jSjl = 1 and {ij) runs over all nearest-neighbor pairs. The 
exchange interaction Jjj takes on either +J or — J with the same probability of 1/2. Using the 
algorithm, we examine the ground-state stiffness of the model using a periodic-antiperiodic 
boundary method and an open-boundary-twist method. 

In section 2, we will give two types of GA for the Heisenberg model. First, we will show that 
the triadic crossover in the Ising model also works well in the Heisenberg model (GAI). Then, 
combining it with a parameter-free genetic algorithm in optimization problems, we will give 
a useful algorithm (GAH). In section 3, using the GAII, the domain- wall energy AE^ of the 
model will be calculated up to a moderate size of the lattice (L < 12) and the stiffness constant 
9 will be given in good accuracy in both the periodic-antiperiodic and open-boundary-twist 
methods. Section 4 will be devoted to conclusions. 

2. Genetic Algorithm 

We develop two GAs in the Heisenberg model: one is the GA proposed by Pal in the Ising 
model ;^ the other is that proposed in other optimization problems. 
2.1 A Triadic Crossover (GAI) 

We first consider Pal's GA, which can be explained briefly as follows. One starts with a 
population of Mj = 2^ random configurations (individuals) of the Ising spins {a^''^ = ±1} (k = 
1, 2, 3, ■ • ■ , Mj), which are linearly arranged in a ring. Then two neighbors from the population 
are taken {parents PI and P2) and two offspring are created using a triadic crossover: For 
each offspring, (i) a different individual {mask or reference R) is selected from the population, 
and the spin of the offspring {(tI^^} is put as af^'' = (t\^^^ or af^'' = cr^^^^ if {(^i^^ crl^^^) = 1 
or -1, respectively; (ii) some fractions fi of the spins of the offspring are reversed {mutation); 
(iii) the offspring is optimized using some local optimization method. Each offspring competes 
with only one parent: the one that is more similar to the offspring. The parent is replaced 
by the offspring if the energy of the parent is not lower than that of the offspring. After this 
step is repeated uoMi times {vq is called a generation period), the population is halved to 
save computer CPU time. From each pair of neighbors, the individual with a higher energy is 
eliminated. This procedure progresses until only four individuals remain in the last stage. One 
merit of this algorithm is that, using the triadic crossover, one can mix the spin configurations 
of PI and P2 over the lattice. Another merit is its ability to slow the loss of population diversity. 
Both are indispensable for the efficient search for the ground state. 

The algorithm might be applied to the Heisenberg model because its components, which 
are characteristic of the Ising model, are the crossover rule and the local optimization method. 
We might replace the equality of the crossover rule with an inequality and the local optimiza- 
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tion method with one appropriate for the Heisenberg model. We start with the population of 
individuals of the Heisenberg spins {s[''^}. We generalize the crossover rule as 

' " \ Sf» for r<"l < , * ' 

where r^-^^^ = {sf^^ sf^^), A; = PI or P27 We use a spin quench method for local optimization; 
all spins are aligned successively in the directions of their local fields. This procedure is 
repeated many times {spin quench step Ng). The similarity of two individuals, k and I, is 
measured with a distance AS^ i between them: 




Actually, ASk,i = when the two individuals are equivalent apart from a uniform rotation, 
although ASk,i ~ ^2/3 (= AS'oo)^ when they are independent. Consequently, all compo- 
nents of the Pal algorithm have been prepared in the Heisenberg model and the algorithm is 
applicable to it. We designated this GA as GAL The CPU time tj of the GAI is estimated as 

ti = AiXi/qX 2(2Mj - 1) X iVg X (lattice size A^), (4) 

where Aj is a constant that is independent of other parameters. Therefore, the CPU time 
depends on three parameters, Mj, z/q, and Ng, in addition to the lattice size N. In other words, 
we must determine the values of those parameters as well as that of the mutation fraction 
fi to optimize the algorithm. Fixing ^ = 0.1, we test the GAI in the itJ Heisenberg model 
on the L X L X L lattice with periodic boundary conditions. For a given bond configuration 
(sample), we apply the GAI 10 times (n = 1,2, 3, • • • ,10) using different initial populations 
and obtain the lowest energy En for each time. We tentatively assume the ground-state energy 
as Eg = mm{En} and presume that we can succeed in getting the ground-state energy when 
{En — Eg)/\Eg\ < 10~^. We made such a test for Ng = 16 samples; results for L < 12 are 
presented in Table I. We see that search efficiency is improved considerably compared with 
that of the conventional spin quench (SQ) method.^ However, to maintain efficiency, we need 
twice or more of the CPU time when L ^ L + 1. One usually encounters this difficulty in the 
optimization problem of complex systems. Here it will be enhanced by eliminating individuals 
at every 'Vq", which degrades the diversity of the population. It is apparent that a large 
period uq slows the loss. In fact, as presented in Table I, search efficiency is improved when 
uq 2uo, but twice the CPU time is necessary. 

In summary, the triadic crossover of the Ising model is applicable to the Heisenberg model; 
furthermore, it improves the ground-state search efficiency considerably. We should optimize 
many parameters for application to different systems. 



3/13 



J. Phys. Soc. Jpn. 



Full Paper 



Table I. Ratios and their standard variances of success in searching for the ground-state energy of 
the ±J model on L X L X L lattice using the GAL where M^, vq, and Nq respectively represent 
the populations in the initial stage, the generation period, and the quenching steps. These are 
estimated for 10 runs for 16 samples. 



L{i^o,N,)\Mi 


32 


64 


128 


256 


512 


1024 


8 (16,200) 


0.34(.23) 


0.52(.23) 


0.74(.23) 


0.90(.20) 






10 (16,200) 




0.59(.31) 


0.75(.30) 


0.87(.19) 






(32,200) 




0.70(.31) 


0.86(.24) 


0.99(.03) 






11 (32,200) 






0.51(.20) 


0.81(.25) 


0.92(.12) 




(64,200) 






0.67(.21) 


0.89(.14) 


0.98(.04) 




(64,400) 






0.87(.15) 


0.97(.04) 


0.99(.03) 




12 (64,400) 








0.81(.23) 


0.95(.09) 


0.98(.04) 



2.2 Parameter- free genetic algorithm (GAII) 

Next, we consider the "Parameter-free genetic algorithm" (PfGA) of Sawai and cowork- 
ers. ^-'^'^ The motivation of the PfGA is to reduce the number of parameters to be determined 
beforehand to apply them to different problems. Although the PfGA reduces such parameters, 
it proves its efficiency in a benchmark test on ICEO.^^ A characteristic point of the PfGA is to 
consider local populations. The numbers of the individuals in those local populations are not 
the same and vary as the population evolves. Another point is the asymmetry of the mutation 
between two offspring, which is inspired by the disparity theory of evolution by Furukawa et 
.12,13 £qj, Qj^g offspring, the mutation fraction /j, is given as a random number, although no 
mutation is applied for the other offspring. The PfGA might also be applied the Heisenberg 
spin system. 

We consider the population S of Np local populations (fc = 1, 2, • • • , Np). The numbers 
of the individuals in those local populations are Uk and the total number of the individuals 
in the population is n(= X^/j^^). As a first step, each local population has two individuals 
nfc = 2 that are constructed randomly and optimized. They evolve as follows: 

(1) Select a local population randomly according to the ratio of Uk/n. 

(2) Select two individuals randomly from S'j. and take out. They consist of a family of S^. That 
is, they are considered as the parents PI and P2 with their energies Epi and Ep2{> Epi). 

(3) Two offspring (children) are produced in the family as follows. 

a. Select one individual, Rl, from different S;(/ 7^ k); generate a child, CI, applying the 
triadic crossover rule described in Sec. 2.1. Generate another child, C2, applying the 
same rule with a different individual, R2, from another Sj,(r 7^ k). 
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b. A mutation algorithm is applied to one cliild (C2)i Clioosc a random number, 

(C'2) 

between (0,1). Attach a different random number, Xj, to each spin S\ . If < Xm, 

(C2) 

S- is replaced by a new one, which is constructed randomly; otherwise, it remains 
unchanged. 

c. The spin quench algorithm is applied for Nq steps to optimize CI and C2. The energies 
of CI and C2 are Eqi and Ec2- If -E-ci > Ec2, then the spin configurations of CI and 
C2 are mutually exchanged. 

(4) Select 1 — 3 members from the family and return to Sj^. The selection rule is as follows. 

a. If Ec2 < Epi, we choose CI, C2 and PI. Here, we consider that PI has superior genes 
and it still has the ability of yielding another good child. 

b. If Ep2 < Eqi, we choose PI. Here, we consider that no superior gene exists in PI or in 
P2. 

c. If Eci < Epi < {Ec2 < or < Ep2 < Ec2), we choose CI and PI. 

d. If Epi < Eci < {Ep2 < Ec2 or £^02 < Ep2), we choose PI and return it to together 
with a new individual. 

In case (a), the number Uk of the individuals of SJ^ increases, whereas in case (b) it 
decreases. We add a new individual to when the number becomes less than two. 

(5) When the best individual {Eci < min{£'p}) is created, we add it to a different local 
population, which is randomly selected. 

(6) We stop searching for the ground state if the best individual is not created during 100 
generations {At = 100). Otherwise, return to (1). 

Figures 1 and 2 show examples of the ground-state search process in a typical sample. In 
Fig. 1, the lowest energies En{t) for eight different populations (eight trials) are presented as 
functions of the generation t. We see that, for each population, stepwise decreases follow an 
abrupt decrease at t 0; its interval At increases with t. Finally, En{t) reaches Eg at t ~ 100 
(further At = 100 generation is necessary to confirm that no more stable state will appear in 
the population). In Fig. 2, the distance ASn,g{t) between the spin configuration {sl^^} with 
En{t) and the ground-state spin configuration is presented in the same search process. Each 
has a value of ASn,g{t) 0.6 (~ 0.7 x AS'oo) at t ~ and changes irregularly with increasing 
t around this value of 0.6 until the system becomes near the ground state. These suggest that 
the system reaches the ground state via local minimum states, the spin configurations of which 
are considerably different from the ground-state spin configuration. That is, the ground-state 
spin configuration will suddenly appear in this search process. 

The PfCA requires only two parameters: the number of the local populations Np and the 
number of quenching steps Ng. In contrast to the GAI, the generation number 1/ results from 
the evolution of the population. It was found that, when Np and Ng are fixed, the average 
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Fig. 1. Lowest energy £'„ for the nth trial in a typical sample of the ±J Heisenberg model as a 
function of the generation t, where Eq is a tentative ground-state energy obtained beforehand, 
applying the same method with a larger number of the local populations of Np ~ 256. Note that 
results for eight trials are shown, one of which (described by *) fails to get the ground-state energy 
because At > 100. 




Fig. 2. Distance in the spin configuration ASn,g as a function of the generation t in the same ground- 
state search process presented in Fig. 1. 



number v slowly increases with the linear size L, e.g., when A'j, = 128 and Nq = 100, u = 
154(16), 149(20), 200(20), 219(26) and 287(30), respectively, for L = 8, 10, 11, 12 and 13. It 
was also found that, for a fixed size L, v decreases slowly as Np and/or Nq increases. Although 
the total number of the individuals in the population S changes as the population evolves and 
often becomes considerably greater than 2Np, the crossover time is 2Np at every generation. 
Therefore, the computer CPU time tn is estimated as 

til = Ml X 2Np X i? X Nq X (lattice size N), (5) 

with i> = 150 ~ 300 for L < 14. The constant An is independent of other parameters and 
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Table II. Ratios and their standard variances of success in searching for the groiind-statc energy of 
the ± J model on L x L x L lattice using the GAII, where Np and Nq respectively denote the local 
population and number of quenching steps. These are estimated for 10 runs for 16 samples. 



HN,)\Np 


16 


32 


64 


128 


256 


512 


1024 


8 (100) 


0.93(.13) 


l.OO(.OO) 












10 (100) 




0.80(.18) 


0.90(.14) 


0.98(.03) 








11 (100) 






0.59(.20) 


0.85(.17) 


0.92(.13) 






(200) 






0.83(.14) 


0.80(.17) 


0.99(.03) 






12 (200) 








0.85(.17) 


0.92(.12) 


l.OO(.OO) 




13 (400) 










0.74(.22) 


0.86(.20) 


0.93(.13) 


14 (400) 












0.71(.19) 


0.81(.23) 



Table III. The ground state energy per spin Ea{L) of the ± J model on the Lx LxL lattice obtained 
by GAII. Eg{= Eg{oo)) is estimated by using data for 8 < X < 13. 



L 


8 


10 11 


12 13 




Ea{L)/J 


-2.0311(2) - 


2.0376(2) -2.0392(3) 


-2.0393(3) -2.0-103(3, 


-2.0132(15) 



All ~ Ai. Wc have performed the same test in the ibJ Heisenberg model; its results are 
presented in Table II. It is readily apparent that tu/ti = 1/4 1/2 for obtaining the ground 
state with the same validity; using GAII, we can extend the treatable lattice size using a 
personal computer from L = 12 to L = 13. We have calculated the ground state energy per 
spin, EciL), for the L x L x L {L < 13) lattice and present them in Table III. Here, the 
parameter set of {Np,Nq) has been chosen such that the search ratio becomes greater than 
0.90. The numbers of samples are about 4000 for smaller lattices and about 500 for larger 
lattices. We have estimated the ground state energy of the model, Eq{= Eg{oo)), by using 
data for 8 < L < 13 with an extrapolation function Eq(L) = Eg + a/L'^, and add it in the 
same table. Note that the value of Eg/ J = —2.0432 ± 0.0015 is a little lower than that of 
Eg/ J ~ —2.0411 estimated by using data for smaller lattices of 4 < L < 11.^ Further studies 
are necessary to settle the value of Eg- 

We have developed a new genetic algorithm (GAII) by combining the triadic crossover of 
the GAI and the parameter free genetic algorithm. The GAII further improves the ground- 
state searching efficiency. It reduces the number of parameters from four (Mj, vq, and ^ in 
the GAI) to two (Np and Nq). This is also the merit of the method because the determination 
of those parameters is an important but tedious task for application to different problems. 
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3. Stiffness of the Ground State 

The most interesting ground-state property of the system is the stiffness of the model at 
T = 0. The stiffness constant 9 is estimated from the domain-wall (DW) energy A.El of the 
L X L X L lattice; 

AEl oc JL^, (6) 

for L ^ oo. When ^ > 0, the SG phase transition occurs at a finite temperature, although no 
phase transition occurs when ^ < 0. A mysterious problem exists, by which 9 depends on the 
estimation method of AEl. Two methods are typically used: One is the periodic-antiperiodic 
(P-AP) method, in which 6 is estimated as 9 = —(0.65 ~ 1.0).^^"^^ However, this value was 
shown to be lattice size range dependent. The other is the open-boundary-twist (OB-Twist) 
method, by which one obtains 9 = 0.5 ~ 1.0 depending on the twisting manner. That 
is, the former method predicts the absence of the phase transition at a finite temperature; 
the latter predicts its presence. Unfortunately, these estimations are given in small lattices of 
L < 8. The question is whether or not the two methods engender the same result for L ^ oo. 
Here, we reexamine 9 in larger lattices using the GAH. 
3.1 P-AP method 

Using this method, one usually considers a simple cubic lattice of L x L x L. Here, we 
consider lattices with different aspect ratios L x L x rL {r > 1) because the value of 9 of 
the system was suggested to be properly estimated in lattices with a large ratio r.^^'^^ We 
consider the DW energy AEj-^l which is defined as the difference in the ground-state energy 
between the two lattices A and B with the same bond distribution but with different boundary 
conditions. That is, for lattice A, a periodic boundary condition is applied for every direction; 
for lattice B, an antiperiodic boundary condition is applied for the z-direction and a periodic 
boundary condition for the x- and y-directions. The GAH is applied to both lattices A and B 
to estimate the respective ground-state energies eI^^ andE^^. TheDW energy is calculated 
as [AEr^L = [eI'^l — E^^^W, where [••■] denotes the sample average. The parameter set of 
{Np,Nq) and the numbers of the samples are the same as described in Sec. 2.2. Figure 3 
shows [A^j-^i] for different r as functions of L in a log-log form, the slope of which gives 
the stiffness constant 9. We see that [AE'^^l] for different r shows different L dependences; 
for r = 1, it decreases with increasing L, but its decrement becomes smaller and seems to 
converge to a finite, nonzero value for L — >^ oo; for r = 2, it is almost independent of L; 
for r > 3, it increases with L and seems to converge to a finite value. That is, the results 
imply that 9 = 0, irrespective of the aspect ratio r. Then using an extrapolation function 
[A£^r,L] = [A£'r,oo] + a/-^ + b/L^ for L > 4, we estimate the convergence value for each of 
r: [AEr,oo]/J = 0.95 ± 0.06,0.49 ± 0.06,0.36 ± 0.06 and 0.22 ± 0.08 for r = 1,2,3 and 4, 
respectively. As expected, results fit well like in a ferromagnetic model: [A^^^^oo] oc [AE^i^coJ/r. 
Therefore, we suggest that = in the P-AP method, contrary to previous estimations of 
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Fig. 3. DW energy AiJ^.L of the L x L x rL lattice using the P-AP method. Lines are fitted ones 
using the function described in the text. 

9 <0. 

3.2 OB-twist method 

We consider the L x L x {rL + 1) lattice with periodic boundary conditions for the x- and 
y-directions and the open boundary condition for the z-direction (the lattice has two opposite 
surfaces, fii and J7rL+i)- We consider a twist energy, which is calculated as follows. We first 
determine the ground state with an energy -Er-.i! then, under the condition that all the spins 
on J7i are fixed, all the spins on f^rL+i are twisted (rotated) at around the z-axis and the 
lowest energy Er^Li4') is calculated; the twist energy is given as the difference between these 
two energies: AEr^L{(j)) = Er^L{<P) — Er,L- Using the GAII with similar conditions to those used 
in the P-AP method, we calculate [AEr^L{4' = for r = 1,2,4. The results are presented 

in Fig. 4 in a log-log form. In contrast to the results in the P-AP method, data for each r seem 
to lie on a straight line, suggesting that [A£'r,L(7r/2)] oc . Their lines' slopes are almost 
identical, which indicates that 9 is independent of the aspect ratio r. Then we fit them using 
a stiffness constant of ^ = 0.62, which was estimated previously in smaller lattices (L < 8) 
with r = 1.^^ The quality of the results is very good, as depicted in Fig. 3. Consequently, we 
expect that 9 ^ 0.62 in the OB-twist method. 

3.3 Remarks related to the stiffness constant 9 

We estimated the stiffness constant 9 in larger lattices using the P-AP and OB-twist 
methods and obtained different values, respectively, of ^ = and 9 ~ 0.62. That is, 9 of this 
model depends on the estimation method, in contrast to a ferromagnetic Heisenberg model 
for which is a universal constant of 9 = d — 2 in the d-dimensional system. We consider 9 in 
those methods. 

In the OB-twist method, the meaning of AEr^L{4>) is clear. We consider the stiffness of 
the ground state of the open system itself; the depth or steepness of the ground-state valley in 
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Fig. 4. Twist energies A£'r,L(7r/2) of the L x L x rL lattice with open boundaries as functions of L. 
Lines are fitted ones using the stiffness exponent 9 = 0.62, which was given using smaller lattices 
{L < 8) with r = 1.17 

the energy landscape of the system. In this case, the twist energy Aii^r,L(0) is surely a lifted 
energy brought by a perturbation around the ground state. That is, the result for 6 ~ 0.62 
suggests that, once the system falls into a ground-state valley at very low temperatures, it 
slightly escapes from the valley. 

On the other hand, in the P-AP method, the calculation of AEr^i comes from an applica- 
tion of the renormalization-group idea; one evaluates the effective coupling Jl(~ JL^) between 
block spins of the linear dimension L generated by renormalization. That is, one assumes that 
Jl ~ AEr^L- However, the meaning of AE^-^l calculated in this method has remained un- 
clear, 20. 23 Q^j, findings of (9 = and [AEr,oo] oc [AEi^^]/r seem to reveal its meaning. We 
consider AE^-^l on the basis of the DW argument in the ferromagnetic Heisenberg model. The 
difference AE^^l arises from the difference in the boundary condition between lattices A and 
B. That is, the difference in the energy between two lattices with different boundary bonds 
{J{i,i),(i.rL)} and {— ^(i,i),(i.rL)}- We first consider the open lattice with { J(i,i),(i,rL) = 0} for 
which the ground-state spin configuration and its energy are described, respectively, as {S^} 
and -E^L- Coming back to lattice some boundary bonds will favor {S^}; others will ob- 
struct it, giving the resulting energy E^j^. We respectively denote the former bonds as right 
(r) -bonds and the latter bonds as wrong(w)-bonds and their numbers as 71^ and riy^. As in the 
ferromagnetic case, this operation will change the energy of the lattice as 0{J /rL) per chain 
along the z-direction for larger L, E^^ — E^ ^ oc —{J/rL) x (n^ — nw). In lattice 5, the roles 
of the r- and tf-bonds are reversed and E^^^ — E^ ^ oc —{J/rL) x (n^, — n^). Then we expect 
AEr^L Of JL^~ /r for L ^ 1, because \nr — ~ L^~ in the d-dimensional lattice. That 
is, = in (i = 3 and [AE'r.oo] [AEi^oc\/r, as was obtained in this study. To examine the 
speculation, we make an additional calculation putting all the surface bonds as r bonds, as in 
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Fig. 5. The DW energy /S.Er^L of the L x L x rL lattice using the Opt- AOpt method. 

the usual DW argument in the ferromagnetic case: i),(i,rL) = 5'«(7n(S^^]^S^?^j^) J}. We call 
this boundary condition an optimized boundary condition and the boundary condition with 
{~'^{i,i),{i,rL)} antioptimized boundary condition, and the calculation of AE'r.L in these 
boundary conditions as an optimized-antioptimized (Opt-AOpt) method. Results of AEr^L 
for r = 1, and 2 are shown in Fig. 5. In fact, AE^^l increases with L for larger L. Interest- 
ingly, their slopes become larger with increasing L, implying that they reach 1 for L ^ oo, 
as suggested by the argument. Further calculations are necessary to estimate the value of 6 
for each r and thereby its r-independent value. Nonetheless the Opt-AOpt method implies 
that > in constrast with ^ = in the P-AP method. We suggest, therefore, that the 6 
in the P-AP method is not a universal constant that describes the effective coupling Ji, but 
a characteristic value that describes the boundary condition. It turns out the results imply 
that a strong correlation exists between the spin configurations on two opposite surfaces of 
the open lattice. 

Hence, we suggest that the result of = in the P-AP method and that oiO ^ 0.62 in the 
OB-twist method indicate the same thing that the ground state of the open system is stable 
against a weak perturbation. 

4. Conclusions 

We developed a genetic algorithm (GA) in the Heisenberg model. We first show that the 
triadic crossover in the Ising model also works well in the Heisenberg model. Then, combining 
it with a parameter-free genetic algorithm, we proposed a useful algorithm for the Heisenberg 
spin-glass model. 

Using the algorithm, we examined the ground-state stiffness of the itJ Heisenberg model 
in three dimensions. Results showed the stiffness constant of = in the periodic-antiperiodic 
method, and that of ~ 0.62 in the open-boundary-twist method. The origin of the difference 
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in 9 between these methods was discussed. We suggested that both results show the same 
thing: the open system's ground state is stable against a weak perturbation. An interesting 
issue whether or not the values of the chirality stiffness constant and the spin glass stiffness 
constant are the same will be discussed in a separate paper. 
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